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Abstract 

I This paper presents original and close to optimal stability conditions linking the time step 

y—{ and the space step, stronger than the CFL criterion: St < C5x°' with a = j^rj : '" a-n integer, for 

f^ some numerical schemes we produce, when solving convection-dominated problems. We test 

CN) this condition numerically and prove that it applies to nonlinear equations under smoothness 

^ assumptions. 

J~z keywords: CFL condition, von Neumann stability, transport equation, Euler equation, Runge- 

Kutta schemes, Adams-Bashforth schemes. 
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I I 1 Introduction 

l^-H In numerical fluid mechanics, many simulations for transport-dominated problems employ explicit 

^L^ second order time discretization schemes, either of Runge-Kutta type [15, 8] or Adams-Bashforth 

(-H [18, 19]. Although widely in use and proved efficient, the stability domains of these order two nu- 

merical schemes (see Fig. 1) exclude the (Oy) axis corresponding to transport problems. Nonethe- 
less, actual experiments [24, 8] show that even in this case, a convergent solution can be obtained. 
If the problem admits a sufficiently smooth, classical solution, the second order time-stepping is 
stable at worst under a condition of type St < C(fa/uijjax)* : where St is the time step, Sx the 
^T) space step, and Umax the maximum velocity of the transport problem. 

^ A close look at the stability condition provided by an analysis of von Neumann type applied to 

^' transport equation provides an explanation. To the best of our knowledge, this result is new -for 

-^ instance, it is not presented in [23] which has collected the state of the art in numerical stability- 

^y-\ despite the fact that it applies to a wide variety of numerical problems. Under some smoothness 

J, • conditions, it readily extends to Burgers equation, incompressible Euler equations, Navier-Stokes 

/-^ equations with a high Reynolds number on domains possibly bounded by walls, and to conservation 

(^ laws. 

^-H For the single step numerical method (i.e. the explicit Euler scheme), a stability result relying 

^ on a similar approach and providing a stability constraint of the type St < C{Sx/u^s.^Y ^^'^ been 

•'^ presented by several authors [12, 17, 21]. The square originates from a completely different kind of 

K^ numerical instability than the usual stability condition for the heat equation with explicit schemes. 

t^ As we will see in this article, it comes from the order of tangency of the stability domain to the 

{Oy) axis and applies only to some first order schemes while for the heat equation it comes from the 

second derivative notwithstanding the order of the scheme. We present the generalization of this 

stability constraint to other schemes. Incidentally, we show that for transport dominated problems 

there exists a direct connection between the order and the stability of numerical schemes. 
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As the numerical viscosity may stabilize the time scheme, this |-CFL-'^ criterion applies essen- 
tially to pseudo-spectral methods and conservative numerical methods [22, 24]. A basic numerical 
experiment allows us to validate our approach. 

The paper is organized as follows: first we recall the definition and the computation of the von 
Neumann stability; then we focus on the linear transport problem, predicting a stability condition of 
the type St < C {5x / uf''^ ' ^'^'^^^> with r an integer, for several schemes; then we construct numerical 
schemes for which such a stability condition appears for r = 1, 2, 3, 4, and corresponds to exponents 
equal to 2, |, | and I; finally we show how this stability criterion extends to nonlinear equations, 
and to multicomponent transport equations (including wave equations). 

2 The von Neumann stability condition 

Let us consider the equation 

dtu^Fu, u(0, O^^uo (2.1) 

where u : M+ x M ^ M, (t, x) h^ u{t, x) and i^ is a linear operator. We denote by ct(^) the symbol 
associated to F, i.e. Fu{t,^) = cr(^)u(i,^) where u n- u stands for the Fourier transform^. 

In the following, we explain how to apply the von Neumann stability analysis as presented in 
[18, 21, 23]. We note Uk ^ u{k6t, •) the approximation at time k6t for k an integer, St denoting the 
time step. We consider we have a spectral discretization, or that all the terms are orthogonally 
reprojected in our discretization space. The scheme can be of Runge-Kutta type, relying on the 
computation of intermediate time steps u/^-j : 



^(0) = Un, U(^i-) 



^ a£i U(i) + (5t ^ 6ft i^ U(j) for 1 < ^ < s', Un+i = u^s') (2.2) 



with {aii)i^i and {bii)i^i well chosen to ensure the accuracy of the integration. 

Or it can be an explicit multi-step (Adams-Bashforth) scheme involving the previous time steps: 

s s 

Un+1 =^^CiUn^i+St^^diFUn-i- (2.3) 

i=0 1=0 

We can also mix these two types of integration schemes: 

£-1 s £-1 s 

"(^) ^ X] '^^^ "(*) + X] "^^^ ""-* + "^^ X ^^* ^ "(*) + ^^ X ^'^^ ^ ""-* for 1 < £ < s', 

i=l 2=0 j=0 i=0 

u„+i=U(s/). (2.4) 

The von Neumann stability analysis consists in isolating a Fourier mode ^ by taking Un{x) = 
(j)n e'^^. Actually, if 6x is the space step, then —f^ < C ^ ij- 

In the case when several previous time samples are necessary, like in the case of an Adams-Bashforth 
scheme, we set 

/ u« \ 



Xn 



(2.5) 



V Un-s / 



Remarking that each time we apply _F to a term in (2.4), we also multiply this term by 5i, it turns 
out that 

X„+i = M{a{^)6t)X^ (2.6) 



-"^CFL stands for the names of the three authors of the founding paper [5]: R. Courant, K. Fricdrichs and H. Lewy 
rn 
defines an isometry on L^ 



^The Fourier transform of a function / £ L-'-(R) is noted /(^) = f_^ /(^) ^ ^^^dx, we recall that / i— > —h= f 
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Figure 1: Von Neumann stability domains for the first four Runge-Kutta and five Adams-Bashforth 
schemes. 



where, setting C, — a{^)5t, M(C) is a (s-t- 1) x (s-l- 1) square matrix whose elements are polynomials 
in C- Note that if F is a differential operator with derivatives of maximal order 7, then |C| < Kj^. 
In the case of hyperbolic equations, 7 is equal to one. 

Let Ao(C), • • ■ , ^s(C) denote the eigenvalues of M{C,). The spectral radius is defined by 



Then 



p(Af(C))=maxo<,<,|A,(C)|. 



p{M{OT < WMicrw < \\M{c)r. 



(2.7) 



(2.8) 



For almost every (, 3^^ > such that Vn > 0, ||Af (C)"|| < i^c PiMiC))" where the constant K(; 
becomes large near the singularities of Af(C). Actually, in numerical experiments, this does not 
play any crucial role (see [21] for a complete discussion on this topic). Hence, overlooking this 
latest point, the von Neumann stability of the scheme (2.4) is assured by: 



Vi,C, |A,(C)|<1 + C<5t 



(2.9) 



with C a positive constant independent of Sx and St. Sometimes, C is taken equal to zero to 
enforce an absolute stability. The assumption (2.9) allows any error Sq to stay bounded after an 
elapsed time T, since: 



llerll = llM(C)^/^*£o|l < i^c (1 + CStp'%o\ 



<Kce^^eo\ 



(2.10) 



The von Neumann stability domain of the scheme (2.4) is given by 5 = {C £ C,p{M{()) < 1}. 
In Fig. 1, 2 and 6, the a;-axis represents the real part of ^ and the y-axis its imaginary part. 
We represent the domain {C G C, |A£(C)| < 1,V^} delimited by the curves {( e C,Xi{C) = e''^,9 G 
[0, 27r[} for < ^ < s. 

On Fig. 1 we plotted such domains for the first four Runge-Kutta schemes and the first five 
Adams-Bashforth schemes. Actually the curves correspond to all the values of C, symbol of the 
operator StF, for which there exists an eigenvalue with modulus equal to one. For order four and 
five Adams-Bashforth schemes, the stability domains only correspond to the semi-disks located on 
the left of the (Oy) axis. The loops on the right of the {Oy) axis do not correspond to any stable 
domain. 



We discretize the differential equation (2.1) with respect to the space variables. We assume 
that u(t,x) = J2k'^kit)ipsx,kix) S Vsx, where Sx is a parameter corresponding to the space step. 
If u = {uk)k, we obtain a discretized version of (2.1) 



dtu = M, 



Sx ^1 



u(0,-) = Uq 



(2.11) 



with, for instance, Mg^ = P^^;^ where P^^; denotes the orthogonal projector on V^x- To ensure the 
stabihty of the simulation, the stability domain with thick line must include the spectrum Sp = 
{\{Msx)} of the matrix M^x- The thick line is due to the term C5t in sup^gg^ p{M{5t£)) < 1 + CSt 
where M from Eq. (2.6) depends on the temporal scheme. 

The behavior of the stability domain along the (Oy) axis indicates how the scheme will be 
stable under the condition p{M{()) < 1 + Cdt -which gives more relevant stability conditions than 
the more classical p{M{Q) < 1- for convection-dominated problems. The next parts of our study 
will be dedicated to finding precise stability conditions on St and 5x in the frame of von Neumann 
stability. 

3 Stability conditions for the transport equation 

The von Neumann stability analysis for the transport equation presents some subtleties which 
explain why Runge-Kutta order two and Adams-Bashforth order two schemes are still used in 
numerical fluid dynamics although the transport operator ia^ for fixed a € M* and ^ G [— y^, ^], 
is located outside the stability domains of these schemes. In this section, we show that what 
matters is the behavior of the stability domain along the (Oy) axis. 

3.1 Accurate theoretical stability condition 

Let us consider the most basic transport equation: 

dtu + adxU^O, with u:R+ xR^R, {t,x)^u{t,x). (3.1) 

Since /'{£,) = i£.f{£,)j the symbol of the operator Fu = —adxU is equal to: a{S,) — —ia£.- 

As explained in the previous section, considering an explicit scheme (2.4), taking Un{x) — (pn e*^^, 

and setting Xn as in equation (2.5), we can write: 

X„+i = AiOXn (3.2) 

with A a matrix whose coefficients arc polynomials in —iaS^St. 

In the case the numerical scheme is of Runge-Kutta type, then A{^) is a polynomial: 

s 

A{0 = Y.M-^^0'St'. (3.3) 

1=0 

The coefficients (/?^) of this polynomial play an important role in our stability analysis. In [18], 
the polynomial g{() — X^fco /^^^^ i'' called the ampUGcation factor. We are able to compute the 
norm of A{^) explicitly: 

1^(01' -E^^'5^''«''^'' (3.4) 

£=0 

with (assuming f3j = for j > s) 

21 

S, = ^(-l)^+^/3,/32,_,. (3.5) 

The von Neumann stability condition |A(^)| < 1 + C5t for all frequencies ^ remaining to the 
computational domain and for a given C, implies that for ^ e I^S' 3a;]' (usually the computational 
domain is rather [— jj:, j^], but we discard tt for simplicity): 

Y,Si5e'a^'e'<^ + '^C5t. (3.6) 

For sake of consistency of the numerical scheme, /?o = /3i = 1 so 5*0 = /^g = 1. Then if 5*1 = • • • = 
Sr-i — and Sr > Q for a given integer r, we can write for small (5i^, 



with (3.6) it implies Sr6t'^''a'^''Sx-^'' < 2CSt, so 5t2r-Q2r^2r ^ q for 5t ^ implies St^ = o(l) i.e. 
as ^ ~ l/(5x, we must have St = o{5x). Hence the equation (3.7) is valid for all the computational 
domain [—j^,j^]- And the stability condition (3.6) is reduced to: 

Sr6t^''a^''6x-^'' < 2C6t (3.8) 



I.e. 

This surprising stability condition is directly linked to the tangency of the stability domain 
S = {C E C,maxi |Ai(C)| < 1} to the vertical axis (Oy). Actually, we have the following theorem: 

Theorem 3.1 (Thick Line Stability Theorem) Consider a numerical time integration of type 
(2.4) with stability domain S bounded near zero by the parameterized curve {C.(0),0 G V(0)} with 
V(0) a real neighborhood of zero. If for some integer r, the Taylor expansion of C, yields: 

C^i{e + o{e)) + T2r9'''' + o{e^^) (3.10) 

with T2r < 0. Then the corresponding stability condition for the transport equation reads: 

Remark that T2r = — ^ where the quantity Sr defined by (3.5) provides (3.9). 

Proof: The amplification factor g{() is given by g{() — Aniax(^(C)) ^ith |Aniax(^(0)| = rnax |Ai(C)| 
and {Ai(C)}o<i<s the eigenvalues of the matrix A{^) from equation (3.2). There is a finite number 
of eigenvalues. Due to the polynomial form of the elements of the matrix A(^), the eigenvalues 
can be written as holomorphic functions: Xi{C) = X]fe>0''^£ C^- Among these eigenvalues Ai(^), 
we consider the one such that |Aio(C)| > |Ai(C)| for all indices i and complex number C in a 
neighborhood of 0. This corresponds to the largest sequence (5o, Si, . . . , Sr, ■ ■ ■) with Si defined 
by (3.5), for the usual order relation on sequences. Then g{() — Xi„{C) is a holomorphic function 
in a neighborhood of 0: 

5(C) = E Z^^^' = 1 + ^ + /^2C' + • • • + M' + ■■■ (3.12) 

where, for consistency reasons, we consider /3o = /3i = 1. 

We already know that for (Se) given by (3.5) satisfying 5*1 = • • • = Sr-i = and Sr > 0, the CFL 

stability condition (3.9) applies. We show that this same condition provides the tangency of the 

stability domain S to the (Oy) axis at zero. 

Near O, let (= p + iq, 

giO = l+p + iq + /32ip + iqf + • • • + Ps{p + iqY + ■■■ (3.13) 

with p and q independent variables close to zero. Then, 

|5(C)|' = {l+P + P2{p'-q^) + ...f + {q + 2P2Pq + ---f- (3.14) 

Looking for the first significant terms of this sum makes 1 and 2p appear. But we do not know 
which is the lowest power of q existing in this sum. Nevertheless, all the terms p^q^ with £,j> 1 
and p with £ > 2 are negligible with respect to p, so we have from (3.13) 

|5(C)|' = \l+p + tq + mqf + ---+M^qr + ...\'+o{p) 

= l + 2p + Srq^'' + o{p) + o{q^'-) (3.15) 



with q the lowest power of q with nonzero coefBcient S^ (given by (3.5)). 

As a result, the curve {|5(C)| = 1} is approximated hy p — —■^q'^^ near the origin. This curve can 
also be parameterized by 6* € [— 7r,7r] in {5(C) — e*^,0 e [— 7r,7r]}. For multistep schemes (see Sec. 
6), it is convenient to express C as a function of 9 and to write it as a Taylor series: 

C-E*^2.-i^''-^ + r2,^''. (3.16) 

Then, for 6* e M close to 

C = t{e + o{9))-^e'^ + o{e'^), (3.17) 

so T2r = ^ 2 • Hence this tangency implies the CFL (3.9). 

D 

Theorem 3.2 An order 2p numerical time integration applied to the transport equation is, at 
worst, stable under the CFL-like condition: 

2p+2 

St<c(-^] (3.18) 

Proof: For an order 2p scheme, we have: 

St"^ St'^P 

Un+l = U„ + StdtUn + —d^Un + ■■■+ -—Ot^U,, + o{dt^P) (3.19) 

2 {2py. 

The transport operator F commutes with dt- So iterating dtu — Fu we obtain d^Un = F^{un). 
Hence equation (3.19) yields the amplification factor: 

3(C) = l + C+? + --- + 7ST+o(C'n (3.20) 



2 (2p) 



with o() gathering the negligible terms under the condition 5t ~ o{5x). In this case, the (/?£) of 
equation (3.12) are given by /3^ = j,. Then for q E [l,p], the coefhcients Sq of the sum (3.4) are 
given by: 

£=0 ^ ^ ' ^ ^' e=a 

Hence, in the worst case regarding the stability, the first nonzero significant term in the sum (3.4) 
is 5p+i(5i2p+2a2p+2^2p+2 ^j^j^ g^^^ ^ g implying the stability condition (3.18). If Sp+i < then a 
linear CFL condition is sufficient. 

3.2 Examples with usual schemes 

We apply our analysis to some popular schemes in fluid dynamics. This provides the following 
stability conditions for some of the most used schemes for transport problem Eq. (3.1): 

• The simplest example is the Euler explicit scheme, order one in time: 

Un+l = Un - St adxUn- (3.22) 

For this scheme, g{Q = 1 + C so r = 1, 5i = 1 and wc find the stability condition: 

5t<2c( — \ . (3.23) 



An improved version of this scheme allows us to construct an order two centered scheme: 



Un+1 = Un - St adxUn^i/2 



(3.24) 



For this scheme, g{() = 1 + ( + ^C^ so r = 2 because 5i = and 5*2 = |. Compared to the 
previous case, the stability is improved: 



St < 2&l^ 



Sx 



4/3 



(3.25) 



• For Runge-Kutta scheme of order 4, we have: 

W„(3) = Un- St a9j;U„(2) 

Un+1 = M„ - fadxUn - faOxU^^i) - f a9a;U„(2) - fada;Un(^3) 



< 



St Q 

Wn(2) — ^n - -2'^ClxUn(l) 



(3.26) 



C 



From the amplification factor g{() = 1 + ^ , ^ 

Si = S2 ~ and 5*3 



C! 



1^ we infer 



1 



-, O4 



(3.27) 



1 

72'"" 576' 

As S3, < 0, our study doesn't apply to this case, and the stability domain, Fig. 1, indicates 
that a classical linear CFL condition has to be satisfied. 



The order 5 Runge-Kutta scheme from [6] page 115 provides the amplification factor g(C) 

1 + C+f + - 



24 "■" 120 



Y^gQ . Therefore it is stable under the condition: 



1/5 



^^ /npoy-^,/,/^ 



6/5 



(T) 



1/5 



4.398 . 



• The order two Adams-Bashforth scheme goes as follows: 

3 1 

Un+1 = Un - -St adxUn + -St ad^Un^l- 



(3.28) 



(3.29) 



So, according to Sec. 2, we consider Xn = 
a pure Fourier mode (/)„ e*^^. We obtain: 



Un 
Un-1 



and we apply the numerical scheme to 



X 



n+l 



1 + IC 



2 




x„ 



(3.30) 



with ( = —iaSt^. 

We compute the eigenvalue of this 2x2 matrix, the characteristic polynomial is given by 
xiY) =y^ — (l + IOy+l. Owing to the fact that St = o{Sx), we have C — > 0. An expansion 
of the larger eigenvalue Yq in terms of powers of ^ provides 



^'o = 1 + C + 



e e c 



o{C 



2 4 

With Q = —ia^, we obtain 

As we want lYol < 1 + CSt, this drives to the following stability condition: 

St < 22/3^1/3 I ''^ ^ 



(3.31) 



(3.32) 



(3.33) 




Figure 2: Von Neumann stability domain for the pseudo- Leap- Frog scheme equation (3.34). 



Therefore, two popular second order schemes, Runge Kutta two (RK2) and Adams-Bashforth 
two (AB2) require a CFL-like condition: St < CSx^^^. The St^^„ is 2^/3 larger for RK2 than for 
AB2, but RK2 necessitates twice more computations than AB2 for each time step. So, regarding 
only the stabihty, AB2 is 2^/^ cheaper than RK2. 

Not all the second order numerical schemes need to satisfy a 4/3-CFL condition. For instance, 
the Leap-Frog scheme calls a usual linear CFL stability condition. The following second order 
scheme is also stable under a linear CFL condition: 



Un+l 



Stadx 



Un +U„-1 



St adxUn 



(3.34) 



Its stability domain is drawn in Fig. 2. The fact that r = 2 with 52 < in Eq. (3.17) is reflected 
by a tangent to {Oy) oriented to the right. 

3.3 Effect of the space discretization 

The space discretization impacts the stability condition (3.18) if it dissipates or creates energy, as do 
the upwind and downwind schemes. Graphically this means that the spectra of these discretizations 
for the transport operator F : u t-^ —adxU are not contained in the {Oy) axis, see Fig. 3. 

In the frame of the von Neumann stability analysis we consider the function u{x) = e'^^, for 
^ £ K. Then, having a closer look at the three academic cases for finite differences, we obtain: 

1. The downwind schemes are always unstable. For the first order downwind scheme (see Fig. 
3 for its spectrum) 



du 
dx 



0{Sx) 



for a < provides the symbol a 



u{x) — u{x — Sx) 
Sx 



P«?^; 



Sx 



(3.35) 



instead of — i a^ in formula (3.3). Combined 



with the Euler scheme for time integration, the amplification factor G{a) 

St f St' 



|G|-l-2a-(l 



Sx 



{l~cos{^Sx)) 



and the error 



et ~ IGl-^'Eo ~ e" 



£o 



1 + Sta becomes 
(3.36) 

(3.37) 



goes unconditionally to +cxi. 



upwind 



downwind 




upwind order 2 
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Figure 3: Spectra of various finite difference scfiemes for space differentiation. The numerical 
method is stable if the spectrum of the discretized differentiation fits into the domain of stability 
of the temporal scheme (as those of Fig. 1). Here, the spectra have been normalized in order to 
have the same vertical size. 



The centered space discretizations satisfy Sp C iM where Sp = {a-(C)iC € [— J^! j^]}- They 
include most of the compact finite difference schemes. For instance the usual centered scheme 



du „ , f 9x "(2; + 5x) — ulx — Sx) ,c^ 

dx^''^'^^ = - — k^ — - = '' 



^i^Sx 



^ — i^Sx 



2Sx 



(3.38) 



has a spectrum given by a 



ijisx) 



which goes along the (Oy) axis, so its distance to the 



domain of stability of the time scheme goes almost the same as in the spectral case. The 
stability results (3.11) presented in this section apply fully to this case with a constant C 
which depends on the space discretization. 

The upwind schemes can be unconditionally unstable if part of their spectra is located on 
the right side of the (Oy) axis (see the spectrum of the order four upwind scheme plotted 
on Fig. 3). If their spectra remain in the left part of the complex plane, then the exponent 
in (3.11) is modified in the following way: assume that the domain of stability of the time 
discretization satisfies 

g{9) =10 + T2g9^^ + o{i9) + o{e^^) (3.39) 

in a neighborhood of 0, with T2q < and q E N; assume that the spectrum of the discretized 
derivative satisfies 

a{9) =i9 + V2p9^P + o{i9) + o{9^p) (3.40) 

with V2p < (i.e. upwind scheme) and p G N. Then the Thick Line Stability condition (3.11) 
becomes: 

• the CFL condition St < CSx, with C a constant independent of St and Sx ii p < q, 

• a mitigation of the nonlinear condition (3.9) 



<;(2p-l) 

St < C(Sa;p<29-i) 



(3.41) 



with C a constant independent of St and Sx, ii p > q. 
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Figure 4: Numerical solution obtained at time T = 1 for three different time steps: 0.97 Sty^nx, 
i^^max and 1.03(5tmax for N = 256 with RK2 (order two Runge-Kutta) with a (5tmax corresponding 
to K = 5 i.e. llurllTy = 5||uo(a;)||Ty- 



The details of the proofs and the numerical tests for these assertions will be presented in a 
further article. Remark that the case p = -l-oo (i.e. switching to a centered finite difference 
scheme) makes the condition (3.11) appear. 

4 Numerical experiment with the Burgers equation 

In order to test our assertions, we proceed to a numerical experiment with the inviscid Burgers 
equation. Although this is a nonlinear equation, we choose initial conditions such that it assimilates 
to a transport equation: the sinusoidal part represents only 1% of the transport amplitude. So, 
technically regarding the stability, it behaves like a transport equation. Then the various Fourier 
modes are naturally activated during the experiment. 



dtu + ud^u = for (t, x) e [0, T] x T, and u(0, •) = uq. 



(4.1) 



In Sec. 5 and 6 we show numerical evidence that stability conditions (3.23), (3.25), (3.33) and 
(3.18) hold for this problem (replacing a by ||u||loo). 

We solve equation (4.1) numerically using a Fourier pseudo-spectral method [2]. The scheme is 
de-aliased by truncation. Most of the time integration methods presented in this paper are tested 
on this classical basic problem. 

The initial condition for the numerical experiment is uo{x) = 10 — 0.1sin(7rx), in a periodic 
domain x € fl = [—1,1]. For t < tmax = IO/tt the equation admits a smooth exact solution 
u{x,t) ~ uo{a), where a = a{x,t) is solution of the equation a — x + uo{a)t = 0. However, the 
numerical solution is only sought for t e [0, 1], in order to satisfy some regularity requirements 
on the solution (see proposition 7.1). To determine the admissibility of the numerical solution, 
we apply a criterion based on the total variation norm (which is expected to be constant): the 
numerical solution u„ has to satisfy ||Mn||Ty < ^||mo(2;)||t\/ with K = 1.1 for all n such that 
nSt<T^ 1. 

The (5tmax we compute, has very little dependence on the divergence criterion K. Actually, 
below (5imax (97%), the numerical solution shows no spurious oscillations, while above it (103%), 
these oscillations create some kind of explosion destroying the profile of the solution completely, 
see Fig. 4. 

The computations are performed for different numbers of grid points, 16 < A'^ < 7758. For each 
N, we find 6tmax by dichotomy with a 0.5% accuracy. The results are represented as 5tmax{N) 
curves in Fig. 5. They evidence the theoretically predicted power law St^ax — CN"" when the 
number of grid points is sufficiently large. The explicit Euler scheme displays a = —2 slope in Log- 
Log scale. The two curves corresponding to the second-order schemes asymptotically both show 
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Figure 5: Maximal time step <5tmax depending on the number of points N for Runge-Kutta schemes 
and Adams-Bashforth schemes, obtained experimentally. 

an asymptotic slope equal to CN^^^ but the constant C is 2^^^ times larger for the Runge-Kutta 
scheme. 

When the order is increasing to 3 and 4 for Runge-Kutta schemes and Adams-Bashforth 
schemes, the slope equals —1. But, while the constant C increases with the order for Runge-Kutta 
(yielding a larger stability domain), it diminishes for Adams-Bashforth schemes with the increasing 
order (see Fig. 1 or e.g., [2]). 

5 Simple 2N-storage numerical schemes with "shrinking CFL" stability 
conditions 

In order to illustrate the phenomenon presented in Sec. 3, we construct numerical schemes having 

stability conditions of the type 5t < (^(x)^"^"^' ^^'^ which only necessitate two time levels to 
be stored in the computer memory. Four of the five schemes presented here need to satisfy this 
stability condition with exponents 2r/(2r— 1) different from 1: 2, |, | and |. All of these numerical 
schemes are of order two, so they show relatively poor consistency given the number of intermediate 
steps. Other efficient low storage schemes can be found in [18] and [11]. 
To solve the equation 

dtu = F{u), (5.1) 
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let us consider the following family of schemes: 

U(i) = Un + apStF{u(o)) 

U{i) == Un + ap_i5tF{u(i,_i)) (5.2) 

Un+l = M„ +ai(5tF(M(p_i)) 

These can also be written: 

Un+i = u„ + ai5tF(un + a2StF{un + a^5tF(un + • • • + ap-i5tF{un + apStF{un)) •••))) (5-3) 
If F is linear, this corresponds to 

Un+l = Un + Pl5tFUn + l32St^F^U„ + fSsSt^F^Un + ■■■+ l3pStPFPu„ (5.4) 

with I3m — Y[T=i '^i- O'wing to F^u — dfu, 

Un+l = Un + PlStdtUn + P25t^dtUn + Mt^dfu^ + ■■■+ ^pStPOfUn- (5.5) 

Here we recognize an expansion similar to the Taylor expansion of the function u„, and we are 
able to tell exactly the order of the scheme for linear equations by comparing the coefficients /3^ 
with those of the Taylor expansion which is provided by: 

Un+l = V -TT^tUn = Wn + StdtUn + -St^d^Un + -Sfdfun + ■■■ (5.6) 

'^ — ' v. lb 

1=0 

and the smallest t such that Pt+i 7^ l/(^+ 1)! indicates the order of the scheme. Remark that this 
holds only if F is linear or if the order of the scheme is less or equal to two. The interest of such 
schemes is that the coefficients ai are easily deduced from the 13 1. 

We assume that i^ is a convection operator. Using the stability analysis Sec. 3, we know that 
the values of 

5, = I3j - 2(it-ifit+i + 2/3,^2^+2 - . . . (5.7) 

provide the stability condition. We verify the validity of this stability condition using the numerical 
test from Sec. 4. 

For a given p in (5.3), maximizing the number of Si equal to zero leads to the following schemes 
of order two -except the first one - and stability conditions: 

• with j3i — \ and /?^ = for £ > 2, this is the Eulcr explicit scheme: 

Un+l = Un + 6tFUn (5.8) 

;32 7^ 5 so it is of order 1, and Si ^ 1 implies 6t < 2C (^)^ 

• with /3i — 1, 132 — 1/2 and f3i ~ ioi £ > 3, this is a second order Runge-Kutta scheme: 

U„+l = Un+ 5tF{Un + -StFUn) (5.9) 

^3 7^ i SO it is of order 2, and 5*2 = 1/4 implies (3.25) St < 2C^I^ i^)'^'^ ■ 

• with /3i = 1, /32 = 1/2, /33 = 1/8 and /3^ = for i? > 4, it is an order two numerical scheme 
(/33 + 1/6), 

(scheme 3) m„+i = m„ + dtF{un + -5tF{un + -StFun)) (5.10) 
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m 


1 


2 


3 


4 


5 


6 


7 


lira 


1 


1 

2 


1 

8 


3-2^2 
8 


5\/5-ll 
64 


26-15\/3 
16 


1 1 7a 
64 ' 128 


ur- 


1 


1.587... 


2.297... 


2.997. . . 


3.687. . . 


3.395... 


5.045. . . 



Table 1: Coefficients /?„ for different m. In the expression for m == 7, a is a real solution of the 
equation a^ — 9a^ — a + 1 = 0. 



and as 5*1 = ^2 = and 5*3 = 1/64, we have the stability condition 



(5.11) 



• the schemes verifying f3i = ioi £ > 5, and 5*1 = 5*2 = ^3 = are given hy f3i — 1, (32 = 1/2, 
k ~ ^ 4 ^ and /34 = ^^g ^ . If wc choose the minus sign for ji^ and /?4, this means: 

(scheme 4) m„+i = u„ + 5tF{un + -6tF{un -\ — StF{un H — StFun))) (5.12) 

It is a second order scheme and has to satisfy the CFL-like stability condition 



6t< 



2CV/WSx^'/' 

PI 



(5.13) 



in the general case, we consider I3q — 1, /3i — 1, /32, • • • , /3m, and for 1 < (, < m — 1, Si — Q. 
As Sm = Pm -^ 0' t^'2 schemes resulting from this system of equations has to satisfy: 



5t< 



2C_ 

w 



6x 



(5.14) 



For Prn positive and minimum, it results the constants indicated in Table 1. 



On the other hand, if we impose the order to be 3 with five nonzero /3^, maximizing the number 
of Si equal to zero provides /3i — \, j32 — 1/2, (i-^ = 1/6, (3^ = 1/24 and (3^ — 1/144. Hence it is 
written 

(scheme 5) u„+i = «„ + StF{un + -7rF{un + —F{un + -ri^(u« + tt -?""«)))) (5.15) 

2 3 4 6 

As Si — Pi — 2l3^l3r, < 0, a classical linear CFL condition 5t < C^ applies. Even, as Pe — l/i\ 
until £ = 4, this scheme is of order 4. 

In Fig. 7, the slopes of stability condition on 5t issued from numerical experiments confirm our 
predictions for these schemes. 

Remark 5.1 Maximizing the tangency of the stability domain to the (Oy) axis is equivalent to 
optimizing the energy conservation scale by scale. This explains why people simulating convection 
dominated problems tend to prefer the Crank-Nicholson scheme Un+i — u„ + ^{Fun + Fu„+i) 
(see [7] for instance) whose stability domain boundary coincides with the (Oy) axis. 
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Figure 6: Von Neumann stability domains for schemes of Runge-Kutta type. 
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Figure 7: Maximal time step ensuring stability for the test case of Sec. 4 with the Runge-Kutta 
schemes stable under the condition St < C5x^^-^ for r = 1,2,3 and 4 whose slops we plot in 
parallel beneath the experimental results. (Ox) axis represents the number of points N — j^ and 
{Oy) the maximal time step (5tinax above which, the numerical solution becomes unstable. 
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6 Adams-Bashforth schemes with "shrinking CFL" stability conditions 

Let us consider an Adams-Bashforth scheme with coefficients (afe): 

K 

Un+l =Un + ^akdtF{Un-k)- (6.f) 

k=0 



The order of scheme (6.1) depends on the sums: 



K 



fc=0 



for < ^^ < i^ 



(6.2) 



(-1)' 



i+ 



J- (sec [18]). Solving the system for 



The scheme has order m, iff for < £ < m ^ 1, Ti 

m = K + 1 provides the Adams-Bashforth scheme of order K + 1 properly speaking. 
The von Neumann stability domain is computed as indicated in Sec. 2. Let 



A"„ — 



Un-l 



and 5tF{ui) — (ui with ^ G 



(6.3) 



V Un^K J 

Then X^ = M{C)X„, with the matrix M(C) G Mk+i{C) given by 



M(C) 



1 + QfoC O^lC 

1 









1 



(6.4) 



The characteristic polynomial is given by 

PiX) = det(X Id - M(C)) = X^+^ -X^ -Y^ akCX 



K 



K-k 



(6.5) 



fe=0 



As the eigenvalues of this polynomial provide the multiplication factor of the scheme (6.1), the 

boundaries of the stability domain are therefore obtained by considering the curve {C G C s.t. 30 G 

[-^,7r],P(e'«) = 0}, i.e. 

e^'^-l 
C=7^7? —^ ee[-T:,^]. (6.6) 






=-ifce ■ 



According to Theorem 3.1, the stability condition depends on the tangency to the imaginary axis 
{Oy) obtained for close to 0. Assuming order one at least (i.e. Tq = 1), a Taylor expansion of 
expression (6.6) provides: 



c-E^ 



E 



(-1)^ 



E (-1)^"^^"" 



p + q = r 
p>l,q>0 



IZ„>inK„=9 



lln>l '*"• 



n 



T 



(6.7) 



The first two elements of this sum are given by 



Ta = -Ti - -, T4 = -Ti - -T2 - T1T2 + -T3 + ^T? + T? 
2 64 62^^ 



(6.8) 



with Ti from (6.2). 

For a given K, maximizing the tangency to (Oy) (i.e. the number m s.t. T21 is equal to for 
I < m) provides the following numerical schemes: 
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Figure 8: Von Neumann stability domains for modified Adams-Bashforth schemes maximizing the 
tangency to the {Oy) axis. 



for K = 1, T2 = impHes ao — ^ and ai = —2, i-c. Adams-Bashforth scheme of order two. 
As T4 = -i, it is stable under the condition (3.33) 6t < 2^/^C^^^ i^)"^^^ ■ 

with three time steps, T2 ^ T4 ^ leads to the scheme we call (ABsch3) with ao = |, 
ai = — I and a2 = \- Given that Ti = —1/2 and T2 = —1/6 7^ 1/3, it is of order two. And 
Tq = —1/12 induces the CFL condition 



6t < 12^/^C^/^ 



Sx 



6/5 



(6.9) 



with four time steps, enforcing T2 = T4 = Tq = yields the scheme (ABsch4) with 

7 21 7 1 
iao,a„a2, as) = {-,--,-,--) 



(6.10) 



As Ti = —1/2 and T2 = 1/10 7^ 1/3, this is also a second order scheme. On the other hand, 
we have Tg 



-^, so this scheme is stable under the condition 



5t < 40^/^(7^/^ 



Sx 



8/7 



(6.11) 



Wc plot the stability domains corresponding to these schemes on Fig. 8, and verify our stability 
predictions with the Burgers equation test Sec. 4. The results of these experiments on Fig. 9 
confirm the predicted stability conditions (3.33), (6.9) and (6.11), but less accurately than for 
Runge-Kutta schemes (3.25), (5.11) and (5.13). 

7 Extension to some nonlinear equations 

We show that these results extend to regular solutions to nonlinear problems such as the incom- 
pressible Euler equations on a domain H. bounded with walls, and scalar conservation laws. We 
proceed in three steps with gradually increasing complexity: 



16 



0,01 



0.001 



0.0001 



le-05 



le-06 



::: ::::::: : :::::: : :: ::::: : 


1 — 










Special A-B scheme 3 — b — [■] 


: ^^^;j;;;^^ ;;- ^ 


Special A 
Adams 


-B scheme 4 — - — : 




; _r^^.;i ;. 














::: ["":["""[""r^i^i^^^"""":""] 












^ ' " " ^"^^ ^'S^" " " " 


...J...... 








..|..|.^.. 


i-;'^^^^^!^^ 


.......... 






: \" 


........ 


- \ \ \ vrW^t''^^ 


^ 












' '';;"; "% 


X 


^w^ 


th. 








;;; ;;;;;;:-ri;;;;r;rr;;;;;T ■"j^:'^ 


!^ 


^ 


:::.:: : y.:\:: 


'.'.y.yy:. 




% 


■"Ot ' 


!"!""! 






■^y>^^>^ : : 


; ; : 




^ 


2S 


m 


-.-. ...... : :::::: : :: ::::: : ::::: : : ::::: 




•fiiLl. 


: : rryTT::': i yviyrrv; r-^T'vr^ -.-:.] 




■■ : nmnl : hhmi -rm- 





10 



100 



1000 



10000 



Figure 9: Maximal time step insuring the stability, obtained experimentally with the 1-D Burgers 
equation (7.15) for Adams-Bashforth schemes. It evidences three slopes: 5tmax — CSx^^^^^^~^'' 
with r = 2, 3 and 4, plus a CFL condition slope, plotted in parallel with dotted lines. (Ox) axis 
represents the number of points N — j^ and {Oy) the maximal time step (5imax above which, the 
numerical solution becomes unstable. 

• First we consider the transport equation with non-constant velocity on bounded domains. 
Hence we step outside the strict frame of von Neumann stability analysis. 

• Then we study the simplest nonlinear equation involving transport: the ID Burgers equation, 
and we show that the previous results still hold true under a smoothness condition. 

• Then we transpose our results to the scalar conservation laws and the incompressible Euler 
equations on a domain D, possibly bounded by walls. 

7.1 Transport by a variable velocity 

The transport of a scalar by a divergence-free velocity u on an open set SI C M'' with regular 
boundaries satisfies the equation: 



dtO + u(x) • Ve* == for X € fl, t € [0, T], 
with div(u) = for X g SI, u(x) • n = for x e dil. 



(7.1) 



In order to generalize the stability analysis to this case, we need the following lemma which 
corresponds to v = [6, ... ,6) and w — {ip, . . . ,(p) m lemma 7.2: 

Lemma 7.1 Let 9,(p : fl ^ M., u : SI — > M"* such that div(u) =0 on SI, and u • n = on dfl then: 

{9, u • V(p)i2(o) = -(u • V0, ip)L^n) (7.2) 

Equivalently, we have: 

{e,M-ve)mn)^Q (7.3) 

The computations using the skew-symmetry relationship leads to the same stability conditions as 
those relying on complex numbers in Sec. 3 under the following assumptions regarding the space 
discretization: 
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Assumption 7.1 The space discretization conserves the skew-symmetry of the equation, i.e. with 
the notations of equation (2.11) 

{Msx0,9) = 0, which is equivalent to (MsxO,^) == -{6,Msx'(p), V^,^ e Vsx- (7.4) 

Such conservative discretizations are presented in sonic computational fluid mechanics publications 
such as [22] for instance. 

Assumption 7.2 The discretization is sufficiently regular to enforce 

\\MsJ\\<C^^, yeeVsx. (7.5) 

This assumption is satisfied with C ~ ||u||loo for almost all the discretizations. One would need 
special properties to avoid this to happen. 

Let F{d) — Psx{u ■ V9) with P^^; the orthogonal projector onto the space of discretization Vsx- 
In the next sections Sec. 7.2, 7.3 and 7.4, even if the operator F is not linear, we omit the projector 
¥sx since it does not change the computations we present because it can be set or removed when 
needed: 

{VsxO,Vsx^) - (PsxO,^) = {e,VsxV), for 9,^ e Vsx- (7.6) 

For the Runge-Kutta scheme (2.2), we find the following expression for 6n+i- 

k 

0„+i=;^ftJfF*(0„) (7.7) 

4=0 

Starting from this expression and according to lemma 7.1 along with assumption 7.1, 

lpi(0\FUf)^\ -/° '^ z + 3 = 2£+l for leN 

(F(6'„),F (0„))i2(o)-| ^_^^,_,||^,(^^^^||2^ if ^ + J = 2^ for leN ' ^^■^> 

We compute the L^ norm of 9n+i as a function of the L^ norm of 6'„. From (7.7), (7.8) and under 
the assumption 7.1, we have: 

fe 

\\en+l\\h=Y.SiSt''\\F'iOn)\\h (7.9) 

£=0 

with (Si) given by (3.5), i.e. 

min(£,fe-£) 
j = -inin(£,fe-£) 

For consistency needs of the numerical scheme, we must have 5o = 1. On the other hand let us 
suppose that ^i = 5*2 = • • • = Sr-i = and Sr > 0. Under the assumption 7.2, for 0„ G Vsx, 

\\F^e„)U.<\\n\\lJ^Ml^, (7.11) 



and knowing that for x > — 1, 
we derive from (7.9): 



yrT^<i + |, (7.12) 



/^Sr\\u\\r^+o{5t)Y'\\e„\\,. < (i+ (^: 



||e„+i|L2 < ( 1 + ^SrMf^ + o{st)) ||e„|L. < ( 1 + ( ^^^^ijuiir^ + 0(1) ) St ] \\e„\\^. 



(7.13) 



where o() gathers all the negligible terms. 

Let us note a = ||u||ioo, then the numerical scheme (2.2) is stable for small perturbations under 
the condition: 

Si<c{^)^\ (7.14) 

Hence the results obtained in the von Neumann stability framework remain valid in the case of the 
convection by a variable velocity on a bounded domain. This is still a linear equation but outside 
the von Neumann stability analysis framework which assumes a periodic or unbounded domain. 

7.2 The Burgers equation 

In order to clarify the role of the nonlinearity and validate our analysis under smoothness conditions 
on the solution, we have a look at the simplest nonlinear case, the one-dimensional inviscid Burgers 
equation: 

dtu + ud^u = Q for (i,a;) e [0,T] x M, u(0, •)=uo. (7.15) 

In order to infer the numerical stability for this problem, we linearize it. Assume m„ is a discretized 
version of the solution u in time and in space. As proposed in [10], we consider a perturbed solution 
Un+Sn- Under regularity assumptions on m, each time discretization will involve a specific evolution 
equation on e„. 

Actually, the small error e„ that we introduce corresponds to oscillations at the smallest scales 
in space Vsx- This instability propagates and may increase at each time step. In the following, we 
demonstrate that under CFL-like conditions similar to those of Sec. 3, the L? norm of the small 
error £„ is amplified in a limited way: 

Ikn+ilU^ <{l + C5t)\\en\\L^ (7.16) 

where C is a constant that neither depends on 5x nor on St. Thus, after an elapsed time T, the 
error increases at most exponentially as a function of the time: 

\\eto+T\\L- < (1 + C5tf"' lletolU^ < e^^lktolU^ (7-17) 

As dtu ~ —udxU, dfu ~ ^^ Xau"^ {dxu)"'^ . . . (9^^^u)"*-i + (— l)^u^9fu, so we remark a kind of 
equivalence between the space regularity and the time regularity. If d^u E L°° then d^u E L°° . In 
the general case, for Runge-Kutta schemes (2.2), we have for < £ < s, 

i-i e-i 



i=0 i=0 



and u„+i + e„+i = U(s) + e^s), 



t-i i-i 

i=0 i=0 



and En+i = £ 



is)- 



Proposition 7.1 Consider a solution u of the Burgers equation (7.15) s-times differentiable such 
that ||9|u||l°°(Rx[o.t]) < +00. Under the condition St = o{6x), a stability error Sn+i in the explicit 
scheme (7.19), small enough at the initial time: HeolU^ = o{5x^i'^) can he expressed as 

s 

e„+i = £„ + ^ /3i (5f < a* e„ + St e„ d^Un + i?„ (7.20) 

i=l 

with ||i?„||i2 = o((5i||e„||i2). The coefficients {j3i) derive from scheme (2.2) similarly as those in 
(3.3). 
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proof: All the terms we have to deal with are projections in the space discretization Vg^- In order 

to simplify the notation, we omit this projection that we assume orthogonal, as for the Galerkin 

methods [24, 1]. 

We prove that e(^) can be put under the form (7.20) by recurrence on i — ... s. 

As £(0) = £n, the assertion is true for £ = 0. 

Let us assume the assertion true for i from to £ — 1 : 

i 

£(i) == £n + ^ I3(i)j 5t^ u^ disn + a(i) St e„ d^u^ + R(i) (7-21) 

with ||i?(j)||^2 = o((5i||£„||i2). The coefficients {P(i)j) correspond to the partial step i of the Runge- 
Kutta scheme distant by au\5t from the time n5t. Remark that au) = 1- Then, given that 

i-i i-i 

^-^ ( ' ... \ 

= £n + ^au \^ fi{i)j5t^u^n9isn + a{i)5tendxUn + R(i) (7.22) 

i=0 \3 = 1 J 

+St ^ bii (U(i)5a;£(j) + £(i)9a;U(i) + £(i)9a;e(i)) • 

Knowing that 

i 
dx£(i) = dxSn + ^ /3(i)jSt^ {dx{ul^)di£n + <9^+^£„) + a{i)5t (dxSndxUn + Sndlun) + dxR(^t), 

(7.23) 
then 

i-i t-i / / * \ 

R(t) = ^ au R(i) + St"^ bu U(i) ^ l3(i)jSt^dx{ui)di£n + St {dxE-nd^Un + £ndlu„) + dxR(i) 

1=0 i = \ \J = 1 / 

+e(i)9i:e(i) + (e{,)9a:M(i) - £ndxU„) + (u(i) - u„) ^ /3(^i)j5t^uidi'^^e„ J . (7.24) 

j=i / 

Now, we need to show that these terms arc o((5i||£„||i2). According to assumption (7.21)and due 
to St — o(Sx), assumption 7.2 provides 

||<^t^a^£„|| < (^^y ||£„|| (7.25) 

Hence £(i) = (1 + o(l))£„ in the sense £(i) = £„ + ?7(i) with ||?y(i)||L2 — odlenlU^) (the stability 
condition being ||£(i)||L2 = (1 + 0{St))\\en\\L^)- As we assumed ||£n||L2 = o{Sx'^/^), then 

\\en\\L^<^^^^o{5x). (7.26) 

As a result, the cross term 5te(i-)dxS(i-) satisfies 

\\5te(^,)dxe(^,)\\L2 < j^\\£{i)\\L^\\£{z)\\L^ ^ o{St\\en\\L^). (7.27) 

As for i < s — 1, 

M(i) ^Un+StB(i'i{Un,dxUn,...,dlUn,St) (7.28) 
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with B a polynomial, ||-B|1l~ is bounded, as well as Wd^BW^^ so |lu(i) — m„||lo°=o(1) and \\dxU{i) — 
dxUn\\L°°=o{l). It allows us to replace ui^\ by u„ in the expansion (7.22), the difference going into 
%), see (7.24). 

Hence, using the fact that eu\,Rfi\ G Vsx the discretization space, ||9^£(i) 11/^2 < '''^ ^ and the 
same for Ruy Let r be an element of the sum Rau then it satisfies: 

Ml^ < ^T(h„|U»,||9,w„|U»,...,||afu„|U~)||£„|U. (7.29) 

with T a polynomial, and p > q + 1. 

Given the fact that 6t = o{6x), we obtain that ||i?(£')||^2 = o((5f||e„||^2). Using the recurrence, we 

obtain the result ioT £ — s i.e. for £„+i. 

Actually, taking into account the orthogonality of Endx^n with £„, we can relax one of the 
assumptions i.e. it is sufficient to have ||eo||L2 — o{Sx), and with the cancellations, it is even only 
necessary that HeolU^ = o((5a;-'-/^). 

D 

Theorem 7.1 If we solve the Burgers equation (7.15) with the numerical scheme (2.2), then for 
a sufficiently regular solution u, the stability condition is provided by: 

where 5t is the time step, Sx the space step, r an integer and Sr a quantity both defined by Eq. 
(3.3), (3.4) and (3.5) and C the constant in the exponential growth of the error: e{t) ^ £06*^* • 

proof: Thanks to proposition 7.1, we are able to write: 

s 
Ikn+lIU^ < ||£„+^ft<5f<a;£„||i2 + St\\endxUn\\L^ + o(5t|| £„ || ^2 ) (7.31) 

1=1 

On the other hand we have ||£n9a;M„||L2 < ||9a;u||L=o ||£„||^2 and since for i,j < s, 

r o(,5t||£„|||,) if ^+J = 2^+l 

{Sfuidle^,dt^uidie^)L-. = { 

[ {-lY-^St^'\\uidisn\\h+o{5t\\er.\\l,) if *+j = 2^ 

(7.32) 
we derive 

s 2s 

||£„+^/3.5f<a;e„||i2 =^5,<5i^||<af£„||i2 +o(5t||£„||i2) (7.33) 

i=i e=o 

with Se given by (3.5). 

Then, as Sq = 1, and ||u^9^£„||i2 < ||u„||f^o IIS"llr.2 



e=o 



6x , 



\e„ + J2l3i6f<dle,,\\l.<Y,sA^) MlUenWh + o{St\\eJl,) (7.34) 



|l£„ + ^/3,<5f<5;£„llL^ < (l + ^E^' 

1=1 \ £=1 



'^(fe) MfL-+o{St)\\\en\\L^ (7.35) 



and finally 



{'^Ip' 



k«+i||L2< {^ + ^y]Si[^^ \\u\\lL + St\\dxu\\L^ + o{St)\ ||£„|U2. (7.36) 
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Let r be the first power in the sum where Sr ^ 0, and let us assume that Sr > 0. Then, the 
stabihty condition ||£n+i||L2 < (1 + CSt) H^kUl^ is reduced to 



1 „ (5i2r-i 



2^'"Sx^r 



llwllfo. < {c - \\d^u\\L^) 



(7.37) 



i.e. the condition (7.30). We recognize the same power law as the one obtained in the linear case 
(3.9). The term — ||9a;u||Loo should usually be discarded since its contribution is external to the 
instability phenomenon and random. 

7.3 Scalar conservation laws 

Scalar conservation laws group equations of the type 



dtu + Y, dx.Mu) = for (x, tjeR'^x [0, T] 
m(0,x) =Uo(x) for xeM'' 



(7.38) 
(7.39) 



with /i : M ^ M diffcrentiable functions and u : M'' ^ M the scalar unknown function. 
A stability analysis of the solution of these equations in the frame of Discontinuous Galerkin 
Rungc-Kutta formulation was presented in [24] for space accuracy of order two and three, with the 
8t < C8x^l^ CFL-like condition, but as the byproduct of a long and rigorous computational process. 
This work was the continuation of [4] where the authors observed that first and second order 
Runge-Kutta methods arc unstable under any linear CFL conditions when the space discretization 
is sufficiently accurate and so does not dissipate too much. In this section, we link their results to 
our analysis and refine the stability criteria. Actually we have the following result: 

Theorem 7.2 Let us apply the numerical scheme (2.2) to solve the equation (7.38). If f G C^"'^ 
andue CP i.e. f(P+'^\u^P'> e L°° and if the (Se) defined by (7.10) satisfy Si = ■ ■ ■ = Sr-'i = and 
Sr > 0, then, given a constant C limiting the exponential growth of the stability error: Et < e'^^egj 
the numerical scheme is conditionally stable under the CFL-like condition: 



5t < 



2C_ 



l/(2r-l) 



5x 



Xtill/;HIU° 



(7.40) 



proof: The proof is more or less the same as for the Burgers case c/part 7.2, using the following 
facts: 



• fi{Un +£n) = fi{Un) + /■ (u«)e„ + o(e„), 

• U(i) -Un = o(l), 

• dxi {f-{un)s) ^ .fi{un)dxi£ for stability analysis, and 

• for all functions rj and £, 



/ d d 
(v,J2J2^^^ (/I("«)/j(u«)9:r,£) 

allowing equalities of the type (7.50). 
Finally, we obtain: 



d 



Y^ fl{un)dx,V, Y fi(^ri)dx,e 



(7.41) 



\i=l 



Wen+iWh ^ [1 + 2Ci6t + o{6t))\\en\\ 



n|lL2 



Sr 5t 



2r 



E 

iG[l,d]'- 



n^K) n^-J^^ 



(7.42) 



L2 
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Then, knowing that for e„ G Vsx, 



iG[l,d]'- \s=l / \s=l / 



L2 



< 



< 



E ii/.K)iu° 

yie[l,d] 



6x^ 



pnlk^ 



3X' 



and neglecting the constant C'l, the von Neumann stabiUty criteria 



kn+i||i2 <(l + 2C5t + o((5t))||e„|| 



L2 



is satisfied if 



i.e. condition (7.40). 



2r 



E ii/K«)iu= 

,iG[l,d] 






(7.43) 



(7.44) 



(7.45) 



7.4 Incompressible Euler equation 

The Euler equations model incompressible fluid flows with no viscous term: 



du 

Ik 



(u-V)u-Vp = 0, divu = 0, for(t,x)e 



X n. 



(7.46) 



The use of the Leray projector P which is the L^-orthogonal projector on the divergence-free space, 
allows us to remove the pressure term: 



du 



+ P[(u- V)u] =0. 



(7.47) 



The stability analysis of this case proceeds somehow as a synthesis of the previous two sections 
Sec. 7.1 and Sec. 7.2. An important property is then the skewness property of the transport term 
(see [9], chapter IV, Lemma 2.1 or [8] for the proof, also used for the stability of the incompressible 
Navier-Stokes equations in [17]): 

Lemma 7.2 Let u, v,w G H^(Vt)'^, H^{U) denoting the Sobolev space on the open set il C M.'^, be 

such that (u • V)v, (u • V)w e L^. If u e fldiv,o(^) = {f e {1^1^)^, div f = on O, f • n = 
on dQ}, then 

(v, (u- V)w)i2(f2) = -((u- V)v,w)i2(s:^). (7.48) 

Corollary 7.1 With the same assumptions as in lemma 7.2, 

(v, (u • V)v)l2(o) = / V • (u • V)v(ix = 0. 



(7.49) 



Considering the scheme (2.2), we introduce a stability error sr£\ at level £. Then, under the 
condition 5t = o{6x) and for e„ small enough, most of the terms appearing in the expression of 
£(£) are negligible with respect to: 

• the terms St' F^(e„) where F„(£„) = P[(u„ • V)£„] and F^ ^ Fn o f „ o ■ ■ ■ o F^ , 

i times 

• the term (5iP[(e„ • V)u„]. 
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Then most of the arguments used in Sec. 7.2 apply with even more accuracy since we have the 
orthogonahty relation: 



{Fn{£n),F^i£n))L2(n) - S ( -1)^-^1 p^ (s\\ 







2 



if i+j ^21+1 for leN 
if i+j = 2£ for ^ e N 



(7.50) 



instead of (7.32). This leads to the following result: 



Proposition 7.2 Assume that the incompressible Euler equations (7.46) have a s-times space- 
differentiable solution u such that || V^uJI^ooQOjTlxn) < +oo, that the discretization conserves the 
skew- symmetry relation (7.48) and thatVe G Vsx, \\Psx Pn{s)\\L^ !i C'l|u||f,oc V k^ ■ Then a stability 
error e small enough at the initial time: HeoHl^ = o{Sx'^''^) remains bounded for t G [0,T] under 
the condition: 

Sr) lllulUo 

with 6t the time step, 6x the space .step, and r and Sr obtained as in (7.10). 



St< 



(7.51) 



This proposition extends to Navier-Stokes equations for high Reynolds number. The incom- 
pressible Navier-Stokes equations are written: 



i9tu + u • Vu - i^Au + Vp = 0, 
divu = 0, a; G 

u(0,a;) = Uo(a::) 



', te[0,T] 



(7.52) 



Using the Leray projector P -the orthogonal projector on divergence-free vector fields- we reduce 
the equation to: 

dtu + P [u • Vu] - i^Au = (7.53) 

Two second order schemes are widely in use for the solution of this equation: the order two 
Runge-Kutta scheme [15, 8] as well as the second order Adams-Bashforth scheme [18, 19]. 



When the Reynolds number Re = 



is sufficiently large, the contribution of the heat 



kernel to the stability vanishes [8] , and the same instability effects as for the incompressible Euler 
equation appear as it was observed in 2D experiments [8]. New tests with boundaries comply the 
stability condition St < (5iniax — C6x '^ for dipole/wall numerical experiments. These results will 
be presented in a forthcoming paper. 

8 Multi-component transport 

We extend the scope of application of the stability conditions (3.9) to other cases with multiple 
derivatives in time, like wave equations, or multiple components, like in some MHD models [7]. 
Let us consider the one dimensional equation: 



dtX = MdxX, with X = 



For example, for X — {u, w)*, and M 



\ Un J 



ue 



and M G A1„(I 



(8.1) 



1 

1 



, we obtain the wave equation dfu — d^u. 
Regarding the general case, we diagonalize the matrix M in C: 

" Ai " 

M = P-^DP, with D = 







An 



•2) 
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Considering Y — PX, the equation dfY — DdjY has physical meaning only if A^ € K for all I. 
Under this form all the components are independent. Therefore all our results on the transport 
equation apply to this case taking a = max^ jAfj. 



When the matrix M cannot be diagonalized, like in the case M 
the second component is independent from the first component: 

dtui = dxUi + dxU2 

dtU2 = dxU2 



1 1 
1 



we remark that 



(8.3) 



then the term dxU2 in the first equation plays the role of a source term. 

In the case when there are several space variables: 

dtX ^ M^d^^X + M2d,,X + ■■■ + M„a,„X (8.4) 

applying a von Neumann stability analysis, we obtain: 

dtX = (Ml iii + M2 i6 + • • • + M„ iin) X. (8.5) 

We consider Af (C) = Mi ^i + M2 6 + • ' ' + ^n ?n, and diagonalize M(^) = P{C)~'^D{(,)P{£,). As 
previously, taking i^ = P{S)X, we obtain the stability constraint (3.11) as in the scalar case. 

9 Conclusion 

The stability CFL-like conditions presented in this paper may be encountered in many simulations 
of convection-dominated problems using explicit numerical schemes. Although based on a classical 
von Neumann stability analysis, this kind of stability analysis is not performed elsewhere. 

Two arguments support our approach. First we explain some "CFL shrinking" effects for second 
order schemes already in use: people remarked that they had to take C — ?• in the usual linear CFL 
condition 5t < CSx. Secondly, we predict some exotic CFL conditions St < C6x^^-^ for certain 
Runge-Kutta schemes and Adams-Bashforth schemes which optimize the energy conservation. 
Numerical tests validate these predictions. 

We showed why increasing the temporal order of a scheme increases the stability. We even 
linked the order of a scheme, its stability and the tangency of its stability domain to the (Oy) 
axis in the von Neumann stability analysis. Nevertheless, the numerical viscosity may erase these 
instability effects especially when using an upwind scheme [4] . 

We extended the domain of application of these results to different equations, including equa- 
tions on bounded domains, non linear equations, and equations with multiple derivatives in time. 
These extensions assume smoothness properties for the solution. This smoothness assumption 
restrains the frame of application to a rather limited area. Nevertheless, this clear exposing of 
actual numerical artifacts plus the correlate accurate stability conditions should be useful to a 
wide community, in particular to those who perform numerical simulations of turbulent flows with 
spectral codes. 
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